
# ------------------------------------------------------------------------------
# Replication Script for:
# "Losing Touch: The Rhetorical Cost of Governing"
# 
# Author: Frederik Hjorth, University of Copenhagen
# 
# This script accompanies the manuscript and reproduces most main analyses,
# figures, and tables as reported in the paper. Remaining analyses are 
# reproduced in the scripts 'topics.R' and 'experiment.R'.
#
# Citation:
# Hjorth, F. (forthcoming). Losing Touch: The Rhetorical Cost of Governing. 
# *Comparative Political Studies*.
#
# Please direct questions or comments to: fh@ifs.ku.dk
#
# Last updated: May 14, 2025
# ------------------------------------------------------------------------------

setwd("~/Dropbox/research/losingtouch/replication") # Remember to set your working directory


library(tidyverse)
library(scales)
library(fixest)
library(modelsummary)
library(viridis)
library(marginaleffects)
library(wbstats)

# Define colors
show_col(viridis(9))
vircol1 <- viridis(9)[1]
vircol2 <- viridis(9)[2]
vircol3 <- viridis(9)[3]
vircol4 <- viridis(9)[4]
vircol5 <- viridis(9)[5]
vircol6 <- viridis(9)[6]
vircol7 <- viridis(9)[7]

# Load speeches data ----
pdf <- read_rds("11-pqmpg.rds")

# samp <- pdf %>% 
#   sample_n(1000) %>% 
#   select(speaker,speechdate,nwords,lix,rarewords_share,ingovt)
# 
# plot(samp$nwords,samp$rarewords_share)

# Preprocess ----
pdf2 <- pdf %>% 
  filter(rarewords_share <= .20 & lix<150) %>% 
  mutate(minister=ifelse(ingovt==FALSE | is.na(ingovt),0,1),
         coarsegovttimer=case_when(is.na(time2firstgovt) ~ "Never",
                                   time2firstgovt<0 ~ "Pre",
                                   minister == TRUE ~ "During",
                                   time2firstgovt>0 & time2lastgovt<0 & minister==FALSE ~ "Break",
                                   time2lastgovt>0 ~ "Post"),
         govtpartyregmp=ifelse(partyingovt==1 & minister==0,1,0)) %>% 
  mutate(rarrev=-100*rarewords_share,
         lixrev=-1*lix,
         rarrevrs=rescale(rarrev,to=c(0,100)),
         lixrevrs=rescale(lixrev,to=c(0,100)))

#how many outliers?
1-nrow(pdf2)/nrow(pdf)

#look at a slice
pdf2 %>% filter(speaker=="Poul Nyrup Rasmussen" & speechdate=="1997-10-21")

pdf2$coarsegovttimer <- fct_relevel(pdf2$coarsegovttimer,"Never","Pre")
table(pdf2$coarsegovttimer)
table(pdf2$minister)


#Find earliest 'break' day for repeat govt. members
firstbreakdaydf <- pdf2 %>% 
  group_by(speaker) %>% 
  filter(coarsegovttimer=="Break") %>% 
  summarise(firstbreakday=min(speechdate))

#merge back in
pdf2 <- pdf2 %>% 
  left_join(firstbreakdaydf,by="speaker")



# Visualize distribution of DVs ----
ggplot(pdf2,aes(x=lixrevrs)) +
  geom_histogram(color="gray80",fill=vircol4,bins = 50) +
  theme_bw() +
  labs(x="LIX measure",y="")
ggsave("figures/histlix.pdf",width=4,height=4)

ggplot(pdf2,aes(x=rarrevrs)) +
  geom_histogram(color="gray80",fill=vircol4,bins = 50) +
  theme_bw() +
  labs(x="Rare words measure",y="")
ggsave("figures/histrar.pdf",width=4,height=4)





# Prepare data for Sun-Abraham estimator for hetereogeneous treatment effects ----

#remove obs from after end of first govt. spell
pdf3 <- pdf2 %>% 
  filter(is.na(firstbreakday) | speechdate<firstbreakday) %>% #takes out speakers from after first spell ends
  filter(coarsegovttimer!="Post") #takes out speakers with only one observed spell post-govt.

nrow(pdf3)/nrow(pdf2) #this takes out 9 pct. of observations

pdf3sunab <- pdf3 %>% 
  mutate(firstgovtyear=ifelse(is.na(firstgovtdate),10000,lubridate::year(firstgovtdate)),
         years2firstgovt=round(time2firstgovt)) %>% 
  select(speaker,speechdate,lixrevrs,rarrevrs,minister,seniority,firstgovtyear,years2firstgovt,year,coarsegovttimer,partyingovt) 

#what is the range of period observations?
table(pdf3sunab$years2firstgovt)
table(pdf3sunab$years2firstgovt)[27]

#how many observations in govt are within first 6 years?
sum(table(pdf3sunab$years2firstgovt)[22:28])/sum(table(pdf3sunab$years2firstgovt)[22:31])
# 93 pct.

#verify that no one is observed during 2nd+ spell
pdf3 %>% filter(minister==TRUE & coarsegovttimer=="Post")
pdf3 %>% filter(minister==FALSE & coarsegovttimer=="Post")

head(pdf3sunab)

#rm(pdf)
#rm(pdf2)
#rm(pdf3)

# Estimate Sun-Abraham diff-in-diff models ----
Sys.time()
res_sunab_lix <- feols(lixrevrs~ sunab(firstgovtyear,year) | speaker + year, cluster = "speaker^year", pdf3sunab)
Sys.time()
Sys.time()
res_sunab_rar <- feols(rarrevrs~ sunab(firstgovtyear,year) | speaker + year, cluster = "speaker^year", pdf3sunab)
Sys.time()

summary(res_sunab_lix)

write_rds(res_sunab_lix,"12-sunab-lix.rds")
write_rds(res_sunab_rar,"12-sunab-rar.rds")
#read back in to save time
res_sunab_lix <- read_rds("12-sunab-lix.rds")
res_sunab_rar <- read_rds("12-sunab-rar.rds")

beepr::beep(14)

(att_lix <- summary(res_sunab_lix,agg="ATT"))
(att_rar <- summary(res_sunab_rar,agg="ATT"))


# print latex table
atttextab <- modelsummary(list("LIX measure"=att_lix,"Rare words measure"=att_rar),output = "latex",
             stars=TRUE,coef_map = c(ATT="ATT of govt. membership"),
             gof_omit = "BIC|AIC",
             title="ATTs for event study models \\label{tab:sunabatts}")

atttextab <-  str_replace_all(atttextab,"& X","& $\\\\checkmark$")
atttextab <-  str_replace_all(atttextab,fixed("speaker\\textasciicircum{}year"),"Speaker-Year")
atttextab <-  str_replace_all(atttextab," speaker"," Speaker")
atttextab <-  str_replace_all(atttextab," year"," Year")

print(atttextab)

write_lines(atttextab,file="tables/sunabatts.tex")

## measure effect sizes in DV sds ----
abs(coef(att_lix)[2]) / sd(pdf3$lixrevrs) 
abs(coef(att_rar)[2]) / sd(pdf3$rarrevrs)


sunabdfl <- broom::tidy(res_sunab_lix)
sunabdfr <- broom::tidy(res_sunab_rar)

sunabdf <- bind_rows(sunabdfl,sunabdfr) %>% 
  filter(term!="seniority") %>% 
  separate_wider_delim(term,names=c("var","lag"),delim="::") %>% 
  mutate(dv=rep(c("LIX measure","Rare words measure"),each=31),
         lag=as.numeric(lag)) %>% 
  filter(lag > -11 & lag < 10) %>% 
  mutate(cilo95=estimate-1.96*std.error,cihi95=estimate+1.96*std.error,
         cilo90=estimate-1.65*std.error,cihi90=estimate+1.65*std.error)

sunabdfref <- tibble(lag=rep(-1,2),estimate=rep(0,2),dv=c("LIX measure","Rare words measure"))

sunabdf <- bind_rows(sunabdf,sunabdfref)

#plot estimates
ggplot(sunabdf,aes(x=lag,y=estimate,group=dv,color=dv)) +
  geom_hline(yintercept = 0,linetype="dashed") +
  geom_vline(xintercept = -0.5,linetype="dotted") +
  geom_errorbar(aes(ymin=cilo95,ymax=cihi95),width=0,position = position_dodge(width=.4),linewidth=.5) +
  geom_errorbar(aes(ymin=cilo90,ymax=cihi90),width=0,position = position_dodge(width=.4),linewidth=1.0) +
  geom_line(position = position_dodge(width=.4)) +
  geom_point(position = position_dodge(width=.4)) +
  scale_color_manual(values=c(vircol3,vircol7)) +
  annotate("text",label="LIX",x=1.8,y=-.9,color=vircol3) +
  annotate("text",label="Rare words",x=1.6,y=-4.1,color=vircol7) +
  labs(x="Time to first year in government",y="Simplicity",color="") +
  theme_bw() +
  theme(legend.position = "none")

ggsave("figures/sunabplot.pdf",width=6,height=4)
ggsave("figures/sunabplot_pres.pdf",width=7,height=4)


# Robustness: add time-varying controls ----

#first get data on time-varying characteristics

#unemployment
dkunemp <- wb_data(indicator="SL.UEM.TOTL.ZS",start_date=min(pdf2$year),end_date=max(pdf2$year),country="DK") %>% 
  transmute(year=as.numeric(date),unemp=as.numeric(SL.UEM.TOTL.ZS)) 
dkcovid0 <- read_csv("99-owid-coviddeaths.csv")
dkcovid1 <- dkcovid0 %>% 
  mutate(year=year(Day)) %>%
  rename(coviddeaths=`Weekly deaths per million people`) %>%
  group_by(year) %>%
  summarise(coviddeaths=mean(coviddeaths,na.rm=TRUE)) 
dkcovid <- tibble(year=1997:2019,coviddeaths=0) %>% 
  bind_rows(dkcovid1)

#merge into pdf3sunab
pdf3sunab <- pdf3sunab %>% 
  left_join(dkunemp,by="year") %>% 
  left_join(dkcovid,by="year")

Sys.time()
res_sunab_lix_timevar <- feols(lixrevrs~seniority + unemp + coviddeaths + sunab(firstgovtyear,year) | speaker, cluster = "speaker^year", pdf3sunab)
Sys.time()
Sys.time()
res_sunab_rar_timevar <- feols(rarrevrs~seniority + unemp + coviddeaths + sunab(firstgovtyear,year) | speaker, cluster = "speaker^year", pdf3sunab)
Sys.time()

(att_lix_timevar <- summary(res_sunab_lix_timevar,agg="ATT"))
(att_rar_timevar <- summary(res_sunab_rar_timevar,agg="ATT"))

beepr::beep(20)

atttextab_timevar <- modelsummary(list("LIX measure"=att_lix_timevar,"Rare words measure"=att_rar_timevar),output = "latex",
                          stars=TRUE,coef_map = c(ATT="ATT of govt. membership",seniority="Seniority",unemp="Unemployment rate",coviddeaths="Covid deaths (weekly per M)"),
                          gof_omit = "BIC|AIC",
                          title="ATTs for event study models including time-varying controls \\label{tab:sunabatts-timevar}")

atttextab_timevar <-  str_replace_all(atttextab_timevar,"& X","& $\\\\checkmark$")
atttextab_timevar <-  str_replace_all(atttextab_timevar,fixed("speaker\\textasciicircum{}year"),"Speaker-Year")
atttextab_timevar <-  str_replace_all(atttextab_timevar," speaker"," Speaker")

write_lines(atttextab_timevar,file="tables/sunabatts-timevar.tex")


# Before and after govt. ----

pdf2_beforeafter <- pdf2 %>% 
  filter(coarsegovttimer!="Never") %>% 
  mutate(coarsegovttimer2=fct_relevel(fct_collapse(coarsegovttimer,Post=c("Break","Post")),"Pre","During"))

table(pdf2_beforeafter$coarsegovttimer2)

m3 <- feols(lixrevrs~i(coarsegovttimer2,ref="Pre")|speaker + year,
            data=pdf2_beforeafter,cluster=c("speaker","year"))
m3_durref <- feols(lixrevrs~i(coarsegovttimer2,ref="During")|speaker + year,
                   data=pdf2_beforeafter,cluster=c("speaker","year"))
m4 <- feols(rarrevrs~i(coarsegovttimer2,ref="Pre")|speaker + year,
            data=pdf2_beforeafter,cluster=c("speaker","year"))

summary(m3)
summary(m3_durref)
summary(m4)

cgm3t <- broom::tidy(m3) %>% 
  filter(term!="seniority") %>% 
  filter(term!="coarsegovttimer::Break") %>% 
  mutate(cat=c("During govt.","After govt."),
         lwr95=estimate-1.96*std.error,
         upr95=estimate+1.96*std.error,
         lwr90=estimate-1.65*std.error,
         upr90=estimate+1.65*std.error) %>% 
  select(cat,estimate,lwr95,upr95,lwr90,upr90) %>% 
  bind_rows(tibble(cat="Before govt.",estimate=0),.)

cgm4t <- broom::tidy(m4) %>% 
  filter(term!="seniority") %>% 
  filter(term!="coarsegovttimer::Break") %>% 
  mutate(cat=c("During govt.","After govt."),
         lwr95=estimate-1.96*std.error,
         upr95=estimate+1.96*std.error,
         lwr90=estimate-1.65*std.error,
         upr90=estimate+1.65*std.error) %>% 
  select(cat,estimate,lwr95,upr95,lwr90,upr90) %>% 
  bind_rows(tibble(cat="Before govt.",estimate=0),.)

cgmt <- bind_rows(cgm3t,cgm4t) %>% 
  mutate(catfac=fct_inorder(cat),
         dv=rep(c("LIX","Rare words"),each=3))

ggplot(cgmt,aes(x=catfac,y=estimate,group=dv,color=dv)) +
  geom_hline(yintercept = 0,linetype="dashed") +
  geom_errorbar(aes(ymin=lwr95,ymax=upr95),width=0,position = position_dodge(width=.4),size=.5) +
  geom_errorbar(aes(ymin=lwr90,ymax=upr90),width=0,position = position_dodge(width=.4),size=1) +
  geom_line(position = position_dodge(width=.4)) +
  geom_point(position = position_dodge(width=.4)) +
  scale_color_manual(values=c(vircol3,vircol7)) +
  annotate("text",label="LIX",x=1.9,y=-.3,color=vircol3) +
  annotate("text",label="Rare words",x=2.8,y=-1.8,color=vircol7) +
  labs(x="",y="Simplicity",color="") +
  theme_bw() +
  theme(legend.position = "none")

ggsave("figures/coarsegovplot_pres.pdf",width=7,height=3)
ggsave("figures/coarsegovplot.pdf",width=5,height=3)


# Compare to govt. backbenchers and to opposition MPs ----

#TWFE for each group
hptwfe_lix_bb <- feols(lixrevrs~minister | speaker+year,data=pdf2,cluster=c("speaker","year"),subset=pdf2$partyingovt==1)
hptwfe_lix_op <- feols(lixrevrs~minister | speaker+year,data=pdf2,cluster=c("speaker","year"),subset=pdf2$partyingovt==0 | !is.na(pdf2$firstgovtdate))
hptwfe_rar_bb <- feols(rarrevrs~minister | speaker+year,data=pdf2,cluster=c("speaker","year"),subset=pdf2$partyingovt==1)
hptwfe_rar_op <- feols(rarrevrs~minister | speaker+year,data=pdf2,cluster=c("speaker","year"),subset=pdf2$partyingovt==0 | !is.na(pdf2$firstgovtdate))

bbopcoefdf <- bind_rows(broom::tidy(hptwfe_lix_bb),
                          broom::tidy(hptwfe_lix_op),
                          broom::tidy(hptwfe_rar_bb),
                          broom::tidy(hptwfe_rar_op)) %>% 
  filter(term!="(Intercept)") %>% 
  mutate(lwr95=estimate-1.96*std.error,
         upr95=estimate+1.96*std.error,
         lwr90=estimate-1.65*std.error,
         upr90=estimate+1.65*std.error) %>% 
  mutate(dv=rep(c("LIX","Rare words"),each=2),
         type=rep(c("Cf. govt. party MPs","Cf. opposition MPs"),2))

ggplot(bbopcoefdf,aes(x=dv,y=estimate,group=type,color=type)) +
  geom_hline(yintercept=0,linetype="dashed") +
  geom_point(position = position_dodge(width=.3)) +
  geom_errorbar(aes(ymin=lwr90,ymax=upr90),position = position_dodge(width=.3),width=0,size=1) +
  geom_errorbar(aes(ymin=lwr95,ymax=upr95),position = position_dodge(width=.3),width=0,size=.6) +
  theme_bw() +
  #  scale_color_viridis_d(begin=.2,end=.75) +
  scale_color_manual(values=c(vircol3,vircol7)) +
  labs(x="",y="Coefficient on govt. membership") +
  theme(legend.position = "none") +
  annotate("text",label="Cf. govt.\nparty MPs",x=.75,y=-1.18,color=vircol3) +
  annotate("text",label="Cf. opposition\n MPs",x=1.12,y=-0.821,color=vircol7,hjust=0) 

ggsave("figures/bbopplot_pres.pdf",width=7,height=3)
ggsave("figures/bbopplot.pdf",width=6,height=4)


# Fit TWFE OLS models ----

m1a <- feols(lixrevrs~minister,data=pdf2,cluster=c("speaker","year"))
m1b <- feols(lixrevrs~minister | speaker+year,data=pdf2,cluster=c("speaker","year"))
m1c <- feols(lixrevrs~minister + seniority | speaker+year,data=pdf2,cluster=c("speaker","year"))

m2a <- feols(rarrevrs~minister,data=pdf2,cluster=c("speaker","year"))
m2b <- feols(rarrevrs~minister | speaker+year,data=pdf2,cluster=c("speaker","year"))
m2c <- feols(rarrevrs~minister + seniority | speaker+year,data=pdf2,cluster=c("speaker","year"))

modelsummary::modelsummary(list(m1a,m1b,m1c,m2a,m2b,m2c))

setFixest_dict(c(lixrevrs="Clarity (LIX measure)",rarrevrs="Clarity (rare words measure)",
                 ministerTRUE="Govt. member",seniority="Seniority",
                 speaker="Speaker",year="Year"))

etable(m1a,m1b,m1c,m2a,m2b,m2c,
       file = "tables/mainregtab.tex",
       digits = 2,digits.stats = 3,
       replace = T,
       title="Regression models of clarity and government membership.",
       label="mainregtab")

# measure effect sizes in DV sds
coef(m1c)[1] / sd(pdf2$lixrevrs)
coef(m2c)[1] / sd(pdf2$rarrevrs)




# Subset to high-profile debates ----

debatedaysdf <- read_rds("11-debatedays.rds")
debatedays <- debatedaysdf$date

#create dummy variable for whether speechdate is in debatedays
pdf2 <- pdf2 %>% 
  mutate(hiprofile=ifelse(speechdate %in% debatedays,1,0))

#TWFE separately for high and low profile debates
hptwfe_lix_hip <- feols(lixrevrs~minister | speaker+year,data=pdf2,cluster=c("speaker","year"),subset=pdf2$hiprofile==1)
hptwfe_lix_lop <- feols(lixrevrs~minister | speaker+year,data=pdf2,cluster=c("speaker","year"),subset=pdf2$hiprofile==0)
hptwfe_rar_hip <- feols(rarrevrs~minister | speaker+year,data=pdf2,cluster=c("speaker","year"),subset=pdf2$hiprofile==1)
hptwfe_rar_lop <- feols(rarrevrs~minister | speaker+year,data=pdf2,cluster=c("speaker","year"),subset=pdf2$hiprofile==0)


debatecoefdf <- bind_rows(broom::tidy(hptwfe_lix_lop),
                          broom::tidy(hptwfe_lix_hip),
                          broom::tidy(hptwfe_rar_lop),
                          broom::tidy(hptwfe_rar_hip)) %>% 
  filter(term!="(Intercept)") %>% 
  mutate(lwr95=estimate-1.96*std.error,
         upr95=estimate+1.96*std.error,
         lwr90=estimate-1.65*std.error,
         upr90=estimate+1.65*std.error) %>% 
  mutate(dv=rep(c("LIX","Rare words"),each=2),
         type=rep(c("All other","High profile debates"),2))

ggplot(debatecoefdf,aes(x=dv,y=estimate,group=type,color=type)) +
  geom_hline(yintercept=0,linetype="dashed") +
  geom_point(position = position_dodge(width=.3)) +
  geom_errorbar(aes(ymin=lwr90,ymax=upr90),position = position_dodge(width=.3),width=0,size=1) +
  geom_errorbar(aes(ymin=lwr95,ymax=upr95),position = position_dodge(width=.3),width=0,size=.6) +
  theme_bw() +
  #  scale_color_viridis_d(begin=.2,end=.75) +
  scale_color_manual(values=c(vircol3,vircol7)) +
  labs(x="",y="Coefficient on govt. membership") +
  theme(legend.position = "none") +
  annotate("text",label="All other",x=.75,y=-.58,color=vircol3) +
  annotate("text",label="High profile debates",x=1.12,y=-1.21,color=vircol7,hjust=0) 

ggsave("figures/debateplot_pres.pdf",width=7,height=3)
ggsave("figures/debateplot.pdf",width=6,height=4)


# Show that estimates are not significantly different
hptwfe_lix_intm <- feols(lixrevrs~minister*hiprofile | speaker+year,data=pdf2,cluster=c("speaker","year"))
hptwfe_rar_intm <- feols(rarrevrs~minister*hiprofile | speaker+year,data=pdf2,cluster=c("speaker","year"))

etable(hptwfe_lix_intm,hptwfe_rar_intm,
       file = "tables/hptwfeintmregtab.tex",
       digits = 2,digits.stats = 3,
       replace = T,
       title="Regression models interacting government membership with debate type.",
       label="hptwfeintmregtab")


# How big a fraction of the days?
length(debatedays)/length(unique(pdf2$speechdate)) # 1.9 pct.

# How big a fraction of the data?
nrow(filter(pdf2,hiprofile==1)) / nrow(pdf2) # 4 pct.


# Separate effects by party -----

ingovtparties <- pdf2 %>% filter(ingovt) %>% pull(party) %>% unique() %>% 
  .[c(1:3,5,7:8)]

#vector of party names translated to English
eningovtnames <- c(
  "Social Democrats",
  "Social Liberal Party",
  "Liberal Party",
  "Conservative People's Party",
  "Socialist People's Party",
  "Liberal Alliance"
)

dk2endf <- tibble(party=ingovtparties,enparty=eningovtnames)


estforparty <- function(partyx){
  cat(str_glue("Getting estimates for {partyx}... \n \n"))
  m1bx <- feols(lixrevrs~minister | speaker+year,data=pdf2,cluster=c("speaker","year"),subset=pdf2$party==partyx)
  m2bx <- feols(rarrevrs~minister | speaker+year,data=pdf2,cluster=c("speaker","year"),subset=pdf2$party==partyx)
  estdfx <- bind_rows(broom::tidy(m1bx),broom::tidy(m2bx)) %>% 
    mutate(dv=rep(c("LIX","Rare words"),1),
           party=partyx)
}

partyspecificests <- map_dfr(ingovtparties,estforparty) %>% 
  arrange(dv,estimate) %>% 
  left_join(dk2endf,by="party")

ggplot(partyspecificests,aes(x=estimate,y=enparty,group=dv,color=dv)) +
  geom_vline(xintercept=0,linetype="dashed",alpha=.5) +
  geom_errorbarh(aes(xmin=estimate-1.96*std.error,xmax=estimate+1.96*std.error),height=0,position=ggstance::position_dodgev(height=.3),linewidth=.5) +
  geom_errorbarh(aes(xmin=estimate-1.65*std.error,xmax=estimate+1.65*std.error),height=0,position=ggstance::position_dodgev(height=.3),linewidth=1) +
  geom_point(position=ggstance::position_dodgev(height=.3)) +
  theme_bw() +
  theme(legend.title = element_blank()) +
  scale_color_manual(values=c(vircol3,vircol7)) +
  labs(x="Estimate",y="")

ggsave("figures/partyestsplot_pres.pdf",width=7,height=3)
ggsave("figures/partyestsplot.pdf",width=6,height=4)


# Plot simplicity over time ----

#calculate overall mean and sd of lixrevrs and rarrevrs
globalmeansds <- pdf2 %>% 
  summarise(meanlix=mean(lixrevrs),sdlix=sd(lixrevrs),
            meanrar=mean(rarrevrs),sdrar=sd(rarrevrs))


avgbyyr <- pdf2 %>% 
  mutate(year=as.numeric(year)) %>% 
  #average simplicity by year
  group_by(year) %>%
  summarise(lixrevrs=mean(lixrevrs,na.rm=TRUE),rarrevrs=mean(rarrevrs,na.rm=TRUE)) %>%
  ungroup()

#line plot for lixrevrs
ggplot(avgbyyr,aes(x=year,y=lixrevrs)) +
  geom_line(color=vircol3) +
  geom_hline(yintercept=globalmeansds$meanlix,linetype="dashed") +
  ylim(c(globalmeansds$meanlix-globalmeansds$sdlix,globalmeansds$meanlix+globalmeansds$sdlix)) +
  theme_bw() +
  labs(x="Year",y="Simplicity (LIX measure)") +
  theme(legend.position = "none")

ggsave("figures/lixtimeplot.pdf",width=8,height=4)

#line plot for rarrevrs
ggplot(avgbyyr,aes(x=year,y=rarrevrs)) +
  geom_line(color=vircol3) +
  geom_hline(yintercept=globalmeansds$meanrar,linetype="dashed") +
  ylim(c(globalmeansds$meanrar-globalmeansds$sdrar,globalmeansds$meanrar+globalmeansds$sdrar)) +
  theme_bw() +
  labs(x="Year",y="Simplicity (rare words measure)") +
  theme(legend.position = "none")

ggsave("figures/rartimeplot.pdf",width=8,height=4)

